Survey of Morphologies Formed in the Wake of an Enslaved Phase-Separation Front 

in Two Dimensions 
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A phase-separation front will leave in its wake a phase- separated morphology that differs markedly 
from homogeneous phase-separation morphologies. For a purely diffusive system such a front, mov- 
ing with constant velocity, will generate very regular, non-equilibrium structures. We present here 
a numerical study of these fronts using a lattice Boltzmann method. In two dimensions these 
structures are regular stripes or droplet arrays. In general the kind and orientation of the selected 
morphology and the size of the domains depends on the speed of the front as well as the composition 
of the material overtaken by the phase-separation front. We present a survey of morphologies as 
a function of these two parameters. We show that the resulting morphologies are initial condition 
dependent. We then examine which of the potential morphologies is the most stable. An analyti- 
cal analysis for symmetrical compositions predicts the transition point from orthogonal to parallel 
stripes. 



I. INTRODUCTION 

Phase-separation is a ubiquitous phenomenon observed 
in a wide variety of systems. The theoretical analysis of 
phase-separation has mostly focused on the case where 
the system is homogeneously quenched, i.e. moved in- 
stantaneously and uniformly from a mixed state to a state 
where the system will separate into different phases. A 
good overview of this theoretical work is provided in the 
book by Onuki[l]. 

In many practically occurring systems phase- 
separation does not occur everywhere at once, but 
rather starts in a specific region and from this place 
successively invades the system. We refer to the surface 
of transition between the mixed and the separated 
regions as the phase-separation front. The resulting 
morphologies formed in the wake of a phase-separation 
front can differ significantly from the structures resulting 
from homogeneous phase-separation. 

Our interest in phase-separation fronts arose from an 
investigation of immersion precipitation membranes [2]. 
In these systems a polymer-solvent mixture, applied 
thinly to a substrate, is immersed in water. As solvent 
leaks into the water and as water enters the polymer- 
solvent mixture phase-separation is induced. It starts 
from the water-polymer /solvent interface until all the 
solvent migrated into the water bath and a porous, asym- 
metric, polymer structure is formed. This structure 
shows a thin initial layer of polymer on the surface that 
will, ideally, show small holes. Below this layer one typ- 
ically finds much larger structures. To understand why 
such structures are formed Akthakul et al. simulated 
immersion precipitation membrane formation in a lat- 
tice Boltzmann framework [2]. Shortly thereafter Zhou 
and Powell examined the same system using a finite dif- 
ference approach [3]. More recently Wang et al. used a 
dissipative particle dynamics method to simulate a the 
effects of varying polymer chain length on the forma- 
tion of immersion precipitation membranes y. However, 
the system proved much too complicated to allow the 



simulations to generate significant insight into the main 
phenomena governing the membrane formation. 

The simulations by Akthakul et al. suggested to us 
the possibility that the main factor controlling the struc- 
ture formation was the dynamics of the phase-separation 
front. However, we found that the dynamics of phase- 
separation fronts in polymer systems was poorly stud- 
ied. This inspired us to investigate the simplest pos- 
sible model for a phase-separation front, i.e. a phase- 
separation front induced by sharp front of a control pa- 
rameter (solvent concentration in the immersion precip- 
itation example) moving with a prescribed speed. 

Phase-separation fronts are also of paramount impor- 
tance in many eutectic alloys, where the physics of the 
front is responsible for the formation and orientation of 
ordered structures, and phase-separation goes hand in 
hand with a solidification problem. Jackson and Hunt 
analyzed the formation of the lamellar (essentially two- 
dimensional) and rod structures sometimes formed in eu- 
tectics. They modeled the dynamics as a steady-state 
diffusion driven phase-separation which occurs directly 
ahead of the solidification front They confirmed the 
earlier observed relation between the solidification front 
speed and the lamellar spacing, and in doing so refined 
some of the earlier theoretical findings. Much of the 
work on eutectic solidification fronts following Jackson 
and Hunt involved adding refinements to their model, 
increasing complexity to make it more like the real al- 
loys under investigation. For instance, the inclusion of a 
convection layer just ahead of the solidification front by 
Verhoeven and Homer [6]. 

Other researchers took an opposite approach; devel- 
oping simple models that could be simulated with nu- 
merical methods. Early work by Ball and Essery simu- 
lating front-induced phase-separation of a binary mix- 
ture noted the remarkable difference between phase- 
separation structures formed by fronts and those formed 
by homogeneous phase separation [7]. Their model con- 
sists of an underlying Ginzburg-Landau free energy, sim- 
ilar to our model in this paper. However, their control 
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parameter is designed to mimic a heat-diffusion process, 
in analogy to the temperature front in eutectic sohdi- 
fication. For slow thermal diffusion they observed the 
lamellar structure familiar to eutecticts researchers, in 
addition to a lamellar structure oriented orthogonally to 
the motion of the temperature front when thermal diffu- 
sion was fast. 

More recently, realistic phase separation fronts induced 
by a control-parameter with its own dynamics, have also 
been studied by Gonnella at alfs*, 'o']. They examined a 
binary fluid where a phase-separation is induced by tem- 
perature change at the walls. The temperature then dif- 
fuses into the finite system, inducing a phase-separation 
front. In this systems the shape and speed of the control 
parameter (in this case the temperature) are space and 
time dependent, complicating the analysis of the result- 
ing structures. 

Another example of a realistic modeling of a phase- 
separation front is given by Kopf et al. who examined 
pattern formation in monolayer transfer for systems with 
substrate mediated condensation [10]. Here the similar 
patterns to the ones predicted in this paper are observed 
in a more complicated system were additional hydrody- 
namic effects lead to a condensation of surfactants in the 
deposition layer. In this case the surfactant concentra- 
tion in the deposited later differs substantially from the 
lipid concentration in the free film so that the concen- 
tration is now dependent on the speed with which the 
phase-separation front advances. 

Hantz and Biro developed a further simplified model 
of a phase-separation front as a moving Gaussian source 
in a diffusive system (llj. By decoupling phase separa- 
tion from a dynamic control parameter, they are able to 
control the phase separation front more directly. In par- 
ticular, Hantz and Biro assembled a rotating front, where 
the front speed is a function of distance from the axis, and 
the direction of motion of the front changes as it sweeps 
across the material. They found, in addition to the ex- 
pected perpendicular and parallel lamella structures at 
the respective slow and fast front speeds, that the lamel- 
lar structures could bend to continue growing perpendic- 
ular to the front. They also observed that the parallel 
lamella formed with width and spacing dependent on the 
speed of the front, and they noted that this was con- 
sistent with experimental observations of Liesegang pat- 
terns. Liesegang patterns are highly ordered structures 
formed in the wake of an electrolyte reaction front in a 
gel[l2|. Antal et al. had earlier used a similar Gaussian 
source front to produce patterns with position and spac- 
ing laws consistent with Liesegang patterns [13]. In an 
earlier paper of ours, we prove analytically that Liesegang 
patterns are reproducible by a similar model [14]. Re- 
fer to the citations contained in that earlier paper for 
other models which have been demonstrated to produce 
Liesegang patterns under the proper conditions. 

An even simpler model of a moving phase separation 
front was used by Furukawa for simulating the formation 
of phase separation induced morphologies [15]. The front 



in his model is an abrupt change in the control parameter 
moving at a constant average speed. Similar to Ball and 
Essery, Furukawa used a Ginzburg-Landau free energy. 
The model by Furukawa is very similar to the one used 
here, with a few minor differences. The implementation 
of our model in numerical simulation, however, is quite 
different as we will show in Sec. IIIII 

In this paper we focus on the simplest possible case: 
we consider a purely diffusive system, as hydrodynamics 
adds additional complexity to this problem. We consider 
fronts moving with a constant speed, since this allows us 
to separate transient phenomena from generic phenom- 
ena, simply by observing the front after a sufficiently long 
time. Furthermore we consider a sharp front, so we do 
not have to consider the effects of the shape of the front. 
For this paper we will focus on two-dimensional systems. 
Despite its simplicity such fronts still exhibit a rich col- 
lection of behaviors, as we will show in this paper. 

This work is an extension of our work on phase sepa- 
ration fronts in one-dimensions. In the one dimensional 
case phase-separation fronts will leave in their wake al- 
ternating domains, and the only remaining question is 
the size of these domains [3, [HI. We were able to show 
that this problem could be solved analytically, at least 
in the limit of small velocities. For two-dimensional sys- 
tems one possible solution are stripes oriented in parallel 
with the front, which are essentially the same as the one- 
dimensional systems observed earlier. 



II. MODEL 

The model we use is essentially the same as the one 
we presented in our earlier paper [16] extended to two di- 
mensions. We therefore give the most important aspects 
in short here: 

To construct a model for phase separation fronts in 
two dimensions we consider a mixture of two materials, 
an ^-type and a S-type, in an incompressible mixture 
such that the total density p = p^(r,t) + psi^^t) is a 
constant. In this paper the position vector r = (x, y) is 
two-dimensional. The order parameter for this system is 
the relative concentration: 



(1) 



For simplicity of the model, we choose a 0^-type mixing 
free energy [17]: 



dx 



c(r,t)0(r,t) + ^M(V<^(r,t))2 



(2) 



The c term, linear in the order parameter, adds a con- 
stant to the chemical potential for spatially homogeneous 
systems. However, in the equation of motion only gra- 
dients of the chemical potential enter the dynamics, so 
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that a constant added to the chemical potential does not 
alter the dynamics of the order parameter. It is included 
here, however, since we will consider different values of c 
across the front, which does influence the dynamics. 

A phase separation front constitutes a spatio-temporal 
change in the control parameter a(r, t) such that the free 
energy at a given location transitions from a mixing state 
(a > 0) with a single minimum, to a separating state 
(a < 0) with two minima. Again for simplicity, we choose 
the transition to be an abrupt jump from a single positive 
mixing value au > to a single negative separating value 
as < 0. The transition moves with constant velocity u, 
and is flat perpendicular to u. The other control param- 
eters are similarly two-valued, with an abrupt transition 
at the front: 

a(r, t) = as ^ {au - a^)© [(r + 1*0 + wt) • u] , 
6(r, t) = hs^ {hu - hs)Q [(r + ro + ut) • u] , 
c(r, t) = cs^ {cm - cs)Q [(r + tq + wt) ■ u] , 
/^(r, t) = K.s^ {i^M - i^s)Q [(r + 1*0 + ut) • u] . (3) 

The mixing and separating values are denoted by sub- 
scripts M and S respectively, and O is the Heaviside 
step function. 

For our diffusive system the equation of motion is 



S,(/)(r,t)=V-[m(r,t)V/i(r,t)] , 



(4) 



where m is the diffusive mobility and /i is the chemical 
potential. The chemical potential is derived from the free 
energy: 

/i(r, t) = ^= t)(/)(r, t) + 6(r, t)(/)(r, tf 

+c(r,t)-/^(r,t)V^^(r,t). (5) 

The equation of motion only considers gradients of the 
full chemical potential. Since fi is itself continuous [16[ 
Fig. 1(b)], we need not be concerned with the compu- 
tational messiness of delta functions which result from 
gradients of the Heaviside function present in the pa- 
rameters of Eq. (j3j). For this model the diffusive mobility 
can take different values across the front: 

m(r, t) = ms + {tum - ms)0 [(r + tq + ut) • u] . (6) 

For the remainder of this paper, space and time depen- 
dence will not be written explicitly, except where needed 
to avoid ambiguity. 

Since we intend to simulate our system with a numer- 
ical method, we will only be able to examine a finite 
system. The way the system is described above, this 
would limit our analysis to times t = Ix/u^ where Ix is 
the length of the simulation in the direction of travel 
of the front. This was a limitation of an earlier, simi- 
lar model by Furukawa[l5j. In turn, this would make it 
costly to investigate the system for large times. To ef- 
fectively look at large times we employ a transformation 
here (as we did in our earlier work ^]), where we keep 



the position of the front stationary in our simulation do- 
main and move the sample with a constant speed u. This 
transformation changes the diffusive equation of motion 
Eq. dH to a drift-diffusion equation of motion: 

dt^{r, t) = V • [-(/)(r, t)u + m(r) V/i(r, t)] . (7) 

Mathematically these two approaches are entirely equiva- 
lent, but for simulation purposes the latter approach has 
the great advantage of allowing us to examine the front 
for long times. 

We rewrite this model in a dimensionless form by con- 
sidering the length, time, and concentration scales of 
spinodal decomposition for a symmetrical system (^^^ = 
0) [Tg*]. The spinodal length is the wavelength of the con- 
centration fluctuation that is most unstable, i.e. phase- 
separating the most rapidly after homogeneous phase 
separation is induced by the quenching a material into 
the spinodal region of the phase diagram. The spinodal 
time is the inverse of the growth rate of those spinodal 
wavelength fluctuations. These are, respectively: 



Xsp = 27r 
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1 tsp 



rusag 



-as 



bs 



(8) 



One of the benefits of this non-dimensionalization is a 
reduction in the number of free parameters of this model 
to the seven following non-dimensional quantities: 



aM ^ 

5^ = 1—5^ 

as OS 



Cm - cs 

as^t^eq 



M = , K = , = -— , U = u—^ . (9) 

ms f^S (Peq A.s 



^sp 



The non-dimensional equation of motion then be- 
comes: 



Ot^ + Vr • (<I>U) 
1 



(10) 



27r2 



where R = r/A^p and T = t/tsp are the discrete non- 
dimensionahzed length and time coordinates. The capi- 
tal script letters are spatially dependent functions of the 
non-dimensional parameters: 



{A,B,C,K,,M} = 



{-1,1,0,1,1}, 
{A,B,C,K,M}, X>Xf 



X<Xf 



.(11) 



The parameters A, 5, M, K are chosen as A = M = 
K = 1 and B = 0. The choice of 5 = is uncon- 
ventional. Its justification is as follows: Since (j) ~ (j)in 
in the mixing region, we can re-expand the free energy 
around (j) = (j)in, retaining only terms to second order. 
This reduces the number of free parameters. As in our 
previous paper [16j, we will restrict ourselves here to the 
case where jusi^^Peq) = l~^M{(t^in)i so that there are no 
long-range diffusion dynamics ahead of the front. Re- 
laxing this restriction will alter the effective at the 
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front and the details of the domain switching. Our choice 
corresponds to C = —A^in. 

This leaves as the only remaining free parameters: the 
initial concentration of the mixed material ^in^ and the 
speed of the advancing front U. Details of the effect of 
changing some of the other non-dimensional parameters 
can be found in our previous work[16]. 



III. SIMULATION METHOD 

To simulate the drift diffusion equation (O we use a lat- 
tice Boltzmann (LB) method, mostly because we intend 
to extend our analysis to hydrodynamic systems where 
lattice Boltzmann methods have been shown to perform 
particularly well. Lattice Boltzmann uses a discretized 
form of the Boltzmann transport equation [18]: 



fi{r + Yi,t+l)- fi{r,t) 



1 



T(r,t) 



[f°-Mr,t)] . (12) 



Time advances in discrete steps {At = 1 is implied 
above), and space is divided into regular cells which tile 
the simulation space. The density distribution functions 
for individual particles are replaced by a discrete set of 
distribution functions /i(r, t) that represent the density 
of particles at position r and time t moving with velocity 
v^. The velocity vectors are restricted such that from any 
given lattice site r, the transformation r r + v^, for 
every index i, always results in a r which lies on a lattice 
site or a boundary site. Following C programming lan- 
guage conventions, the lattice sites are numbered from 
to Ix — 1 in the x direction. 

The zero-order velocity moment of the non-equilibrium 
distribution functions is the order parameter: 



,t). 



(13) 



The choice of the equilibrium moment distributions de- 
termines the equation of motion to be simulated. For a 
drift-diffusion equation, we chose: 

i 

^ /^(r, t)Via = sua(l){r, t) , 

i 

f^{r, t)ViaVi(3 = s/i(r, t) + s'^UaU(3(l){r, t)5al3 . (14) 



E 



The subscripts a and /3 are indices for the spatial dimen- 
sions X and y; for instance Ua represents the magnitude of 
the vector u in the a direction. The time-scaling param- 
eter s introduced here could easily be absorbed into the 
other parameters, but it provides us with a convenient 
dial to select the fastest stable simulation parameters. 

In this transformed reference frame, the control pa- 
rameter front that induces phase-separation is station- 
ary, and the material is advected across the front. We 



choose to align the advection velocity with the x-axis: 
u = {ux^O). The control parameter front is then imple- 
mented by setting the parameters for x < Xf to their 
separating values, and setting the parameters for x > xj 
to their mixing values. For example, the model parame- 
ter a(r,t) from Eq. (|3]) implies the simulation parameter: 



a(r) 



as ^ X < Xf 

dM , X > Xf 



(15) 



Derivation of the equation of motion from these mo- 
ments can be accomplished by a Taylor expansion of the 
lattice Boltzmann equation ([12]) with repeated substitu- 
tion of the unknown functions fi with first order approx- 
imations in terms of the known equilibrium f^. This 
is shown in more detail in our previous paper for the 
one-dimensional system[16], and the analysis in higher 
dimensions remains essentially the same. In terms of the 
simulation parameters, the equation of motion becomes: 



1 



- V • (u(^) + O {d^) = V • (r - 1/2) V/i . (16) 



So we have to identify the mobility m with m = r — 1/2. 
It turns out that this method does not show numerical 
artifacts for rapidly changing mobilities. This suitibility 
for simulating abrupt changes in r was previously used 
to simulate abrupt changes in dielectric properties of a 
medium[19]. 

To fulfill Galilean invariance and isotropy requirements 
on a square lattice, implied by our choice of moment dis- 
tributions Eqs. (p!4|) for an arbitrary advection velocity 
u(r,t)-as would be the case were this simulation cou- 
pled to a hydrodynamic flow-would require the use of 
nine (or seven on a hexagonal lattice) velocity vectors in 
two dimensions. However, because u = i^^^ is fixed and 
aligned to the x-axis we require only five velocities in two 
dimensions: 



(17) 



a so-called D2Q5 LB implementation. For this velocity 
set we use the equilibrium distributions: 



fo = (1 



2 2\ 

uts ) 



2/23 , 



/3 =/4°=W2• 



(18) 

In lattice Boltzmann units the material velocity is —UxS, 
and can be made arbitrarily small by adjusting s. 

The equilibrium distributions contain the chemical po- 
tential /i given by Eq. (j5j). The chemical potential con- 
tains a Laplacian which is evaluated on the discrete lat- 
tice as 



(19) 



5 



where the weights Wij 
trix 




are the elements of the stencil ma- 



(20) 



where the row and columns are numbered { — 1,0,1}. 
This choice of discrete Laplacian with (x± 1, 1) terms 
is less susceptible to certain instabilities than those which 
have only (x ± 1, and (x, y ± 1) terms, allowing us to 
use an almost twice as large an effective time step|20[ at 
very small computational cost. 

We still need to define the boundary conditions. The 
^/-dimension boundary conditions are periodic. The in- 
fiow boundary condition dit x = Ix is straightforward 
with homogeneous material advancing with a constant 
prescribed fiux: 



fi{x = lx-l,y,t^l) = f2{x = lx-l,y,t)^Ua:S(l)i, 



(21) 



To calculate the Laplacian on this boundary we simply 
set (t)(lx,y) = (t>in' 

The outfiow boundary condition is somewhat more 
complicated, since we now have phase-separated material 
that is advected out and we need information about (j) and 
ji from lattice sites that are not part of the simulation 
space. We define our outfiow boundary condition with 
the understanding that it should be neutral wetting to 
all concentration values, should not introduce gradients 
on the chemical potential, and any effect on the morphol- 
ogy should be short ranged compared to the simulation 
size. This is accomplished by simply bouncing back the 
outfiow density distribution after subtracting the mate- 
rial advected out of the simulation: 

f2{x = 0,?/,t + 1) = = 0,?/,t) - UxS(l){x = 0,y,t) . 

(22) 

For the purpose of calculating the Laplacian at the out- 
fiow boundary we set off-lattice concentrations to be the 
same as the boundary (j){—l,y) = (j){0,y). The effect of 
the outfiow boundary condition on the bulk of the sim- 
ulation can be ascertained by observing a normally sta- 
tionary morphology as it is advected out of the simula- 
tion. One choice is a single large circular region of S-type 
material suspended in a bulk which is otherwise entirely 
^-type. Tests such as these (not shown) verify that this 
boundary condition has a very small effect, and is accept- 
able for this model. Further evidence can be found by 



comparing phase-separation structures in Fig. |l(a) and 
Fig. |2(b)| which are nearly identical, as discussed later. 

This fully describes our method which allows us to sim- 
ulate the dynamics of structure generation by phase sep- 
aration fronts. As a final note: while simulations will 
show material moving past a stationary front, in our dis- 
cussion of these simulation results we will always refer 
to the front as moving and the bulk of the material as 
stationary. 




(a)T = 82, 
UT = 10.5 




(b)T = 164, 
UT = 20.9 




(c)T = 966, 
UT = 123.6 



FIG. 1. (Color online) Simulation showing stripe formation in 
the wake of a phase separation front. For this simulation we 
choose a front speed U — 0.128 and a symmetrically mixed 
$ = initial material concentration. The starting location of 
the front is marked as a vertical stripe. The rightmost image 
is the final structure observed after the front has moved far 
into the material. A sample of the final morphology, noted 
by the circled region, will be compared to samples from other 
simulations performed with different parameter values. 



IV. FIRST SURVEY 

Using our lattice Boltzmann method we now investi- 
gate how the phase separation front will infiuence the 
formation of structure. To do this we set up a medium 
size simulation of x = 512 hy y = 1024 lattice points. 
We put the position of the front at x/ = 384 (in lat- 
tice units). The region from x = Otox = Xf under- 
goes homogeneous phase-separation. We use the initial 
condition <l>(r) = + 0.01^, where ^ is random noise 
uniformly distributed in the range [—1, 1], as is typically 
done for homogeneous phase-separation simulations. The 
other simulation parameters are: qm = 1, cis = — 1, 
bM = 0, bs = 1, Cm = -0m, cs = 0, uim = ms = 1/2, 
and hZM = i^s = 2, consistent with the choice of non- 
dimensional parameters made at the end of Sec. [Ill The 
time scaling parameter s = 0.026 is chosen as a max- 
imum numerically stable effective time step for these 
parameters [21]. 

Such a simulation for U = 0.128 is shown in Fig. [TJ 
The current position of the front at any time is eas- 
ily visible as the transition between the black-and-white 
phas e-sep arated region and the gray mixed region. In 
Fig. 1(a) the front has moved only a short distance of 



X = 132 = 10.5 Asp (in units of the spinodal wavelength). 
The position where the front started is marked by a ver- 
tical stripe. The area to the left of the initial location 
of the front has undergone normal spinodal decomposi- 
tion generating a phase-separation morphology typical 
of homogeneous phase-separation. As the front moves 
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(a)T = 33, 
UT = 8.4 




(b)T = 82, 
UT = 20.9 



(c)T = 491, 
UT = 125.7 



FIG. 2. (Color online) Similar to Fig.[Tl this sequence shows a 
simulation resulting in stripe formation by a phase separation 
front moving with speed U — 0.256 into mixed material of 
initial concentration $ = 0. At this faster speed the stripes 
formed by the front are oriented parallel to the front. 



on, new material phase-separates, but we see immedi- 
ately that the structure of these newly formed domains 
is quite different from the domains formed through ho- 
mogeneous phase-separation. In the region between the 
initial front position and the current front position we 
observe a different kind of morphology: the domains 
are oriented orthogonally to the front, and show a va- 
riety of widths. In Fig. |l(b) the front has advanced 



a distance x = 263 = 20.9 A^p. The region of homo- 
geneous phase-separation has noticeably coarsened and 
the newly overtaken material continues to phase-separate 
into a striped morphology. The striped structure, how- 
ever, is not homogeneous. Particularly where the stripes 
are thin, defects c an be seen to traverse into the striped 
domains. In Fig. l(c)| the front has traveled a distance 



X = 1553 = 123.6 Asp. Now all the defects have been 
advected out of the system and the stripes are taking on 
a uniform thickness. This is a stationary solution that 
will persist indefinitely. Note that there is no evidence of 
any further coarsening at the position of the front. 

To examine the effect of front speed we now perform 
a simulation for U = 0.256. The results of this simu- 
lation are shown in Fig. [2j This time the morphology 
formation at the front is qualitatively different. Again 
a regular striped morphology is formed, but it is now 
oriented parallel to the fr ont in agreement with previ- 
ous results [nl UH. In Fig. |2(a) the front has traveled a 



short distance of x = 105 = 8.4 A^p. We observe typ- 
ical homogeneous phase-separation morphology behind 
the original front location, however, where the front has 
traversed there are stripes of somewhat regular widths 
oriented roughly parallel to the front. While the stripe 
widths are fairly uniform, there are still a large number 
of bends in these stripes. In Fig. |2(b)| we see that as time 



progresses new stripes form with fewer sharp bends, but 
the stripe widths do not appreciably change. 

Note that the region of homogeneous phase-separation 
corresponds to the homogeneous phase-separation in 
Fig. |2(a) at the same non-dimensional time. The ini- 
tial noise on the order parameter was identical for both 
simulations and closer examination shows that the result- 
ing phase-separation morphologies are nearly identical. 
This shows that there is little interaction between the re- 
gions of striped morphology and the homogeneous region. 
It also shows that neither the outflow boundary condi- 
tion nor the advection speed significantly influence the 
simulations. This indicates that the coarsening of the re- 
gion which separated under homogeneous conditions only 
slightly affects the striped morphology. 

The width of stripes oriented parallel to the front 
can be understood by considering this as a quasi one- 
dimensional system. We analyzed this situation in an 
earlier paper [16] and found that the wavelength of the 
parallel stripes follows from the front speed as 



L{U) 



! In 2 



(23) 



for very slow fronts {U < 0.001) moving into material 
that has vanishingly low diffusive mobility ahead of the 
front (M 0). For faster speeds this relation breaks 
down, and this theoretical prediction is inappropriate for 
the front speed of U = 0.256 in the example simula- 
tion shown in Fig. [2l However, the observed wavelength 
of these quasi-one-dimensional stripes L||2d = 1-39 with 
M = 1 compares favorably with the measured stripe 
wavelength Lid = 1-36 from the M = simulation re- 
sults in our earlier paper on phase separation fronts in 
one-dimensional systems [16]. 

These two qualitatively different morphologies were 
first described by Furukawa[15l. They were later redis- 
covered by Hantz and Biro [11] and appear also to be 
related to structures formed from eutectic mixtures, al- 
though typically a phase-field formalism is used to de- 
scribe these structures. 

The fact that simply changing the velocity leads to 
a change in orientation of the domains raises the ques- 
tion as to where this transition happens, and surprisingly 
there appears to be no numerical value for the speed at 
which this transitions happens in the literature. Also if 
we change the input composition, we will have to obtain 
stripes that have different width depending on the com- 
position. This more systematic investigation constitutes 
the main contribution of this paper. Next we obtain a 
state-diagram from simulations such as the simulations 
as shown above. We pick a sample of the morphologies- 
shown as a circle in Fig. |l(c)| and Fig. |2(c) -and place 
these sample morphologies in a diagram at a position 
corresponding to and U. We then performed sim- 
ulations for a set of different values of and [/, but 
with all other parameters kept constant. We ran each 
of these simulations until the front had moved approxi- 
mately 1500 lattice sites, or about 4 times the distance 



7 




FIG. 3. (Color online) Morphology phase diagram from simulations started with random initial conditions. Examples from 
this map are shown in Fig. [T] [21 and[4l The bar in the upper left corner has length 10 and height 1 in non-dimensional units. 
The lack of an apparent pattern is due to the strong hysteresis of morphologies. Further explanation is given in Fig. |4] and the 
text. 



from outflow boundary at x = to the front position at 
Xf = 384. We again use a circular secti on near the front, 
as indicated in Fig. |l(c) and Fig. |2(c)[ and put them in 
a ^in/U graph at the appropriate position. This survey 
of resulting mophologies is shown in Fig. [3l 

The result of the survey is initially surprising: while 
there is a clear transition from parallel to orthogonal 
stripes for = at around U ^ 0.2, orthogonal do- 
mains appear to be an anomaly only observed for exactly 
symmetric domains. However, there appears to be an- 
other transition between orientations of asymmetric do- 
mains at around a hundredth of this speed at ^ 0.003. 
The state diagram shows one additional boundary be- 
tween regions of stripes and droplets for more off-critical 
mixtures. These droplet structures appear to be pre- 
ferred for larger speeds. For even larger speeds and more 
off-critical mixtures we see only mixed material, visible 
as gray disks. This means that the speed of a phase- 
separation front moving with constant speed into the 



mixture is smaller than the imposed speed of our front, 
and phase-separation is unable to keep up with our front. 
The transition to a free front, to the accuracy of this 
survey, is unaffected by the initial conditions. For the 
parallel stripes the free front speed is given by(22| 



UiU^in) = (34 + 14^7) (3 - 9^lf . (24) 

Closer examination of this state-diagram reveals more 
unexpected results. There are three examples of orthog- 
onal stripes formed in a sea of parallel stripes at (/7, ^in) 
values of (0.045, 0.05), (0.045, 0.15), and (0.01, 0.35). 
These structures appear to break the prevalent trend of 
their neighbors and it is worth while to consider these 
simulations in some more detail. We will focus here on 
the simulation for {U = 0.01, = 0.35). Three snap- 
shots of this simulation are shown in Fig.HJ The homoge- 
neous spinodal morphology are droplets. However at the 
front a depletion zone favoring white material is formed. 
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(b)T = 2619, 
UT = 29.6 




5565, 
62.9 



FIG. 4. (Color online) Simulation showing strong depletion 
effects. The moving front suddenly appears in a mixed ma- 
terial with volume fraction = 0.35 containing small fluc- 
tuations. Shown in (a), the region far behind the front un- 
dergoes spinodal decomposition and droplet growth as in a 
homogeneous quench. Near the front the formation of a de- 
pletion zone induces domains oriented parallel to the front. 
By chance a defect in this domain facilitates a transition (b) 
to the favored orthogonal stripe morphology shown in (c). 
The circle in (c) shows the region sampled for use in Fig. O 
The front speed is U — 0.0113 for this example. 



The reason for this depletion layer is as follows: we 
chose to set the C-parameter in the mixing region such 
that /i(^in) = in the mixing region. After phase- 
separation we have = ±1) = in the separated 
region and there will be no long range chemical potential 
gradient leading to extended diffusion of material into 
or out of the mixing region. Before the initial phase- 
separation in the separating region the order parameter 
is nearly uniformly and ii{^in) ^ 0. Therefore we get 
some diffusion into and out of the mixing region leading, 
in this example, to a depletion of black material near the 
front. Once the phase-separation is complete, there is no 
longer a difference in the chemical potenial far away from 
the front. 

After the initial phase-separation with the creation of 
the depletion zone we then observe the nucleatio n of a 
first parallel black domain, as can be seen in Fig. 4(a)[ 
So far this scenario is generic, but this simulation is spe- 
cial in that the formation of the first parallel stripe is not 
perfect, but a single defect was created. This defect now 
has a profound effect on the further evolution of the mor- 
phology formation. The next black domain that is formed 
has two defects with an interspersed drop. The follow- 
ing generation of parallel stripes has three drops, but the 
middle drop now maintains contact with the front, form- 
ing the first orthogonal stripe. The defect invades the 
formation of parallel domains leading to a wedge of or- 
thogonal stripes that replace the parallel stripes. After a 
sufficiently long time we are left with a purely orthogonal 



morphology. 

This suggests that instead of our state-diagram as 
shown in Fig. 2] we should associate the orientation of 
the domains with a probability, since the selection is ap- 
parently probabilitstic, depending on the details of the 
homogeneous phase-separation, which in turn depends 
on the initial noise. However, since the appearance of a 
single defect can be sufficient to switch the orientation 
(as shown in Fig. |4]) we expect the probability of find- 
ing a certain morphology to also depend on the system 
size, since it is more likely to develop such a defect in 
a larger system. In the limit of a macroscopic system, 
we would expect that the probability of finding a defect 
would increase significantly, so that it becomes more in- 
teresting to examine which morphology is the preferred 
morphology. 



V. SECOND SURVEY 

To find out what is the preferred morphology we choose 
an initial condition which contains stripes of both paral- 
lel and orthogonal orientation and is consistent with the 
overall composition ^^n- Parallel stripes are generated 
through a nucleation process and such a stripe selects 
its preferred length. For the formation of orthogonal 
strips no nucleation occurs and it is more difficult for 
such stripes to select a preferred length scale. Here we 
design our initial condition with a range of stripe widths 
to allow for easy selection of the preferred width. We 



show this initial condition in Fig. 5(a) 

We construct the initial conditions by selecting a 
square region of the simulation with side length / that 
is half the simulation height. Assuming this region has 
an origin in the lower left, the concentration in the region 
is given by: 



^(x, y) - signum 



sm 



al — {a-\- b)y 



sm 



where a ■ 



Xsp and b ■ 



(25) 

4 Xsp are, respectively, the small- 
est and largest stripe wavelengths initially generated. 
The region described is the lower-center of Fig. 5(a)[ The 
upper-center region is constructed similarly, although 
with a 7r/2 rotation to favor formation of parallel oriented 
stripes. To avoid the possibility of the outflow boundary 
condition interfering with stripe selection dynamics, we 
disconnect the initially phase separated regions from the 
boundary with a region of mixed material at the initial 
volume fraction containing small fluctuations similar to 
the region ahead of the front. 

In Fig. [5(b)] we see that only one of these initial stripes 
is selected to form the first orthogonal stripe. For this 
example the orthogonal stripes are again the preferred 
morphology (Fig. |5(b)] ) and the orthogonal stripes even- 
tually replace the parallel stripes (Fig. |5(c)' ). Depending 
on the volume fraction and the front speed either orthog- 
onal stripes (as shown here), parallel stripes, or droplets 
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(a)r = O^UT = 0.0 




(c)T = 2292, UT = 51.9 




(b)r = 573, UT = 12.9 



(d)T = 9739, UT = 220.4 



FIG. 5. (Color online) Example simulation showing selection 
of favored morphology for an imposed front speed U = 0.0226 
and initial concentration = 0.35. These parameters have 
a favored morphology of stationary stripes, but is close in pa- 
rameter space to the droplet morphology region. The initial 
phase separation configuration (a) allows for the selection of 
the favored morphology without the strong depletion effects 
observed in Fig. 4(a) which sometimes occur after spinodal de- 
composition near a front. The stability of the periodic stripes 
is evident in (c). The final configuration in (d) shows the sta- 
ble stationary stripe morphology, and a circle which outlines 
the region sampled for use in Fig. [6l 



may turn out to be the preferred morphologies. 

With this new initial condition we may now calculate 
a new state diagram. The simulations used in this sur- 
vey use the same non-dimensional parameters as the first 
survey. Here the simulation size has been changed to 
X = 768 hy y = 1024; widened to accommodate the new 
initial conditions. To compensate for the increased com- 
putational cost of a larger simulation we chose to use a 
smaller value for the interfacial free energy cost = 1, 
which increases the effective simulation speed by allow- 
ing a larger effective time step s = 0.079 [20]. This new 
state diagram is shown in Fig. [6l As expected we find 
that for the preferred morphologies there is now a much 
larger range of parameter values for which we find orthog- 
onal stripes. For symmetric mixtures the transition be- 
tween orthogonal and parallel stripes appears unchanged 
up to the accuracy of this survey. As we change the 
from symmetrical {^in = 0) to asymmetrical {^in 0) 
compositions we now observe a continuous transition in 
morphologies. This is because we start with a phase- 



separated morphology in the separating region, and we 
no longer form a depletion layer. Thus parallel or or- 
thogonal stripes are not a priori favored. We see that 
increasing the volume fraction still lowers the speed for 
which we see a transition between orthogonal and parallel 
stripes, albeit in a much less drastic fashion. 

We also find a larger parameter range for which we 
find droplet arrays, particularly for larger speeds and 
more symmetrical volume fractions. We can now write 
a schematic state diagram for the preferred morpholo- 
gies. This is shown in Fig. [71 The region labeled "free 
front (periodic)" is where phase separation lags behind 
the control parameter front. The "free front (single do- 
main)" regions are where initially undifferentiated ma- 
terial will not spontaneously demix without fluctuations 
to induce nucleation. Without fluctuations, instead of 
new domain formation, any existing domains of phase- 
separated material will slowly grow into the mixed ma- 
terial. How structures which may form in these regions 
are effected by the passage of a control parameter front 
is beyond the scope of this paper. The regions under 
the solid curve are based on observations of the preferred 
morphology evidenced in Figs. [3] and [6l Apart from the 
free front region we also show the regions where we ob- 
serve parallel stripes, orthogonal stripes, and droplet ar- 
rays. The former two we have covered in previous sec- 
tions, but the later requires more discussion. The droplet 
structure is observed to initially form near the front with 
little long-range order. As the front progresses, the po- 
sition of the newly forming droplets is influenced by the 
depletion of material caused by the formation of the pre- 
vious drops. The larger the drop formed, the more it de- 
pletes the surrounding region, and the further away the 
next droplets will form. This mechanism causes reorder- 
ing and elimination of small droplets in favor of larger 
droplets. A nucleation condition imposes a maximum 
droplet size for a given front speed and mixed material 
concentration. The result is a droplet structure which 
converges towards a highly ordered array of homogeneous 
size droplets. Whether the droplet array has a preferred 
orientation, or if it is periodic in the reference frame of 
the front, are open questions. 



VI. ANALYSIS OF RESULTS FOR CRITICAL 
CONCENTRATIONS 

We have observed in the surveys that a phase sepa- 
ration front moving into material that is at the critical 
concentration = will form a striped morphology 
that is oriented either parallel or orthogonal to the front. 
The parallel stripes are an essentially one-dimensional 
morphology, and have been discussed in our previous 
paper 0. As mentioned in Sec. |IVl the key property 
of parallel stripes is that for a given set of parameters 
there is a unique stripe wavelength, and the wavelength 
scales L|| oc 1/VU. By contrast, orthogonally oriented 
stripes may form with a wide range of wavelengths for a 
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FIG. 6. (Color online) Morphology phase diagram from simulations started with initial conditions containing a variety of 
morphologies. The initial morphologies include long stripes oriented both parallel and orthogonal to the front with stripe 
wavelengths that range from one to ten spinodal wavelengths, and a region of mixed material with small fluctuations near the 
outflow boundary. An example of the initial condition and final morphology selection is shown in Fig. [5] 



given parameter set. In this section we analytically de- 
termine the stable wavelengths for orthogonal stripes as a 
function of the enslaved front speed. We accomplish this 
in two stages. We first determine the maximum wave- 
length of a stripe before a new stripe nucleates in its cen- 
ter, splitting it. We next find that there is a met ast able 
minimum stripe width, below which coarsening dynam- 
ics can result in stripes which disappear into the front. 
The result is a region of stability for the orthogonal stripe 
morphology, outside of which we predict the favored mor- 
phology is parallel stripes or droplets. 



A. Maximum Orthogonal Stripe Size 

The mechanism responsible for limiting the maximum 
stripe wavelength is the nucleation of an opposite-type 
stripe at its center. This can occur, even without fluctu- 
ations, due to the build-up of a nucleation kernel ahead 



of the enslaved front, similar to what happens in the one- 
dimensional case [16]. 

Formation of stable orthogonal stripes by a moving 
front results in a chemical potential profile that is sta- 
tionary in the reference frame of the font. We analyti- 
cally determine this profile directly from the stationary 
solution {dr^ 0) to the equation of motion: 

= ^V^M - Vr • ($U) . (26) 

The morphology under consideration is an array of highly 
ordered orthogonal stripes of wavelength L± with an A- 
type stripe centered on the X-axis and the phase separa- 
tion front on the F-axis. The front velocity U = ([/, 0) is 
entirely in the X-direction and constant, simplifying the 
gradient term to a derivative of the concentration ^ in 
the X-direction which is highly dominated by the con- 
centration gradient at the front. We approximated this 
term as an abrupt step which we expect to become exact 
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FIG. 7. Morphology state diagram of front induced phase 
separation generated structures in two dimensions. See text 
at the end of Sec. |V]for a detailed description of the regions 
and boundaries. 



in the limit of large stripes and slow fronts. We focus 
on critical mixtures here, since additional complications 
occur for off-critical situations, as mentioned below. For 
critical mixtures, the gradient term then simplifies to 



Vr • (U$) 



dX 



U5{X) Sqr(r/ij 



(27) 



where Sqr(x) = signum[cos(27rx)] is the square-wave 
function. This is the only approximation made in our 
calculation. 

An alternate form of the square- wave function is as an 
infinite series of delta function convolutions. To make the 
series converge faster, each term of the series is centered 
on a crest and extends to include half of each adjacent 
trough: 



Sqix)= J2 



/ n+1/4 n+1/2 

2 J S{x — x) dx — J S{x ■ 

\n-l/4 n-1/2 



x) dx 



Combined with the delta function property 6{x/a) = 
\a\5{x) and the n-dimensional vector delta function iden- 
tity ^n(r) = S{xi)S{x2) ■ • ■S{xn)^ this alternate form al- 
lows Eq. (|26]) to be rewritten 



)dY 



oo / ^' 

V^/i = 27i'^UL^ ^ ^ / ^2(i^±?ey - R) ( 

n+1/2 

- j ^{Ljrey -Yi)dY I , (29) 

n-l/2 



where ey is the F-axis unit vector. The boundary condi- 
tion for this morphology is a vanishing chemical potential 
far ahead and behind the front: ±oo) = 0. The 

Poisson equation V^/(r, tq) = (^2(1*0 — r) has the fun- 
damental solution /(r, ro) = ^ In |ro — r|. This can be 
directly applied to the above equation to arrive at an in- 
finite series integral solution for the chemical potential 
profile: 



/ n+1/4 

^ / r I ~ ^ 

/i(L^R) = itUL^ 2 / In y + (F - dY 

^=-^\n-l/4 

n+1/2 

- In \[x^ -^{y - yf dY I , (30) 



(28) 



n-1/2 



It is easily seen that, by symmetry, this solution satisfies 
the chemical potential boundary condition. The inte- 
grals are straightforward, however the result is somewhat 
lengthy: 



/i(L^R) = nUL^ ^ i (n - F) 



. f {n-Y)/2 
2 arctanh —7^ ; ---7^ — 



2 arccot 



X 



1/4 -F 



-2 arccot 



X 



1/4 -F 



arctanh 



arccot 



n-Y 



X2 + (n-F)2 + l/4 
X 



1/2 -Y 



- arccot 



X 



1 / [X^ + (n + 1/4 - Yf] [X^ + (n - 1/4 - Yf] 
'2 I 16 [X2 + (n + 1/2 - F)2] [X2 + (n - 1/2 - Yf] 



1/4 -F 
(31) 



I 

We verify the analytical chemical potential in Eq. (|3T]) by comparison to the chemical potential measured from 
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FIG. 8. (Color online) Comparison of the analytical chemical 
potential profile of Eq. (|3ip (line) with LBM simulation results 
(circles) for large stripes L± = 70 formed by a slow U = 
0.00235 front. The left figure shows /a along the center line 
of a stripe of B-type material. The right figure shows /x at 
the position of the front ahead of two adjacent stripes. The 
profile in each figure intersect at X/L± — Y/L±_ — 0. 



a LBM simulation of a stripe pair with a large stripe 
wavelength = 70 and slow front speed U = 0.00235. 
The result of this comparison is shown in Fig. [8l and 
demonstrates excellent agreement for large stripes and 
slow fronts. The assumption of Eq. (|27|) is invalid for 
fast fronts, and the X/L± profile becomes asymmetric 
(depressed on the leading edge, and bulging on the trail- 
ing edge) and both profiles show an overall supression 
of the chemical potential from the analytical prediction. 
These additional profiles are not shown here. 

However, to predict the values of L and U where nu- 
cleation occurs we only need the extremum value of the 
chemical potential. This occurs in the center of a stripe 
at the front and is given by 



/i(R = 0) = ttUL^ 



1 /16n2 

2 Vl6^ 



1 



(32) 



+ nln 



32n^ - 6n - 1 
32n^ - 6n + 1 



-2UL^C = /iP^^^ . 



Here C = 0.9159 .. . is a constant known as Catalan's 
constant. Nucleation of a new stripe will occur when the 
chemical potential peak reaches the nucleation chemical 
potential //P^^^ = /i"^^^. 

Due to the symmetry at the center of the stripe, the nu- 
cleation chemical potential for a two-dimensional stripe is 
the same as the one-dimensional switching condition. In 
a previous paper [16] we presented an analytical expres- 
sion for the switching concentration for the special case 
M = where there is negligible diffusive mobility ahead 
of the front. In that case an analytical expression for 
the switching chemical potential can be found. However, 
a nonzero M induces earlier switching due to the pres- 



ence of a nucleation kernel ahead of the front [IGj. With 
the same method we used to verify the switching condi- 
tion in our previous paper, we have measured the one- 
dimensional switching chemical potential /i"^^^ ^ 0.24 
for the non-dimensional parameters used in this two- 
dimensional system. Thus we find that the maximum 
orthogonal stripe wavelength as a function of front speed 
is: 



^max 



0.24 0.131 



2CU 



U 



(33) 



It is interesting to note that the analytical composi- 
tion profile ahead of a solidification front for lamellar 
morphologies of eutectic mixtures presented by Jackson 
and Hunt is also a solution to the full chemical poten- 
tial profile, although in a very different form[5, Eqn. (3)]. 
By using the substitutions C ^ ji^ = Sf3 ^ ^±/4, 
Co = ^ 1/2, V U, d l/27r^ x Y , and 
z IXI, their solution can be written: 



/i(L^R) = 2UL_ 



CO 

E 



1 sin f^) cos(2n7rr)e-2»-l^l , 
V 2 / 

n=i 

(34) 

which is, at least in a preliminary numerical evalua- 
tion, equivalent to our solution, although we have not 
yet shown this analytically. For off-critical situations a 
complication occurs: the inflow material induces a non- 
neutral wetting condition for the orthogonal stripes. This 
can be clearly seen in Fig. IH This situation requires ad- 
ditional considerations that will be discussed elsewhere. 



B. Minimum Orthogonal Stripe Size 

The minimum stripe width is limited by the width of 
stripes which, if a defect occurs, coarsen more quickly 
towards the front than the front moves away. A simple 
qualitative argument for the defect speed can be obtained 
from the dynamical scaling laws [23]. This law describes 
the time evolution of a morphology by stating that at 
later times the structure is statistically similar to that of 
earlier times when scaled with the typical length scale 



(35) 



The constant C is expected to depend on the parameters 
of the system and the details of the morphology. From 
Eq. (|35]) we can define a coarsening speed which provides 
an estimate of how quickly the end of a finger of wave- 
length Lc will recede: 



Uc 



1 



dT 



(36) 



If a finger protrudes from the front into the phase sep- 
arating region, it will coarsen in the same direction as 
the front. If the coarsening speed is faster than the front 
Uc > U ^ the finger will eventually coarsen away com- 
pletely into the mixed material domain. According to 
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FIG. 9. (Color online) Measurement of the dynamical scaling 
constant C in Eq. ()35p using LBM simulations. Simulation 
size is 4096^ lattice sites except for hi — 2.0 which is 1024^. 
The numerical method is similar to that which is outlined in 
Sec, mil for Ux — 0, except with fully periodic boundaries and 
no front. Characteristic length scale was measured by divid- 
ing the system area by the length of the $ = interface, then 
non-dimensionalized by multiplication with a scaling factor 
(2.27/Asp) such that Lc = 1 is the length scale of spinodal 
decomposition. 



Eq. ([36|) smaller fingers coarsen at a greater speed, and 
this sets a minimum wavelength for orthogonal stripes 
generated by a moving front: 

>ic(c/) = 2y^. (37) 

Unfortunately, numerical values of C for different situa- 
tions are hard to come by in the literature. However, it is 
easy to obtain C by measurement during phase ordering 
of a homogeneous quench. We did this for a symmet- 
ric system (i.e. <I>(T = 0) ^0) and obtained a value of 
C ~ 0.555; for details please refer to Fig. [9l 

We expect this coarsening speed to be on the order 
of, but not exactly equal to, the speed with which a sin- 
gle finger coarsens. To obtain a better estimate of the 
coarsening speed of a finger we can measure the speed in 
a simulation. To do that we set up a proto-stripe sim- 
ilar to what is shown in Fig. [TOl We then measure the 
position of the tip of the finger as the first zero-crossing 
of (j) at the original y position of the center of the fin- 
ger. We then vary the front speed slowly (once every 
Alx iterations) to stabilize the position of the finger tip. 
Once the velocity has reached a stationary state (deter- 
mined by the tip speed and the average of the last 100 
tip speed measurements being less than /7/10^) we find 
the coarsening speed of the finger. The results of those 
measurements are shown in Fig. [11] and, as expected, the 
stripe speed is well approximated hy U = 4C^/3L^, with 
C = 0.29. This is approximately a factor of 2 smaller 
than the dynamical scaling constant found in Fig. [9l 




FIG. 10. (Color online) An example of a proto-stripe sta- 
bilized by regulation of front speed. For the color version 
of this plot the range of values are shaded - in the ad- 
ditive RGB color model - red=(/i - flmin) / (l^max - /J^min), 

green=((/) - 0min)/(0macc -0min), bluc^rcd X grecu. A legend 
is shown at the right of the figure. Equi-chemical potential 
lines are superimposed. The front speed is increased when 
the stripe shrinks and decreased when the stripe grows until 
a stationary system is found; see text for methodology. Initial 
conditions are similar to the final state (shown), with an ini- 
tial front speed U — (0.2L^)^ although the final U was found 
to be insensitive to the initial value. The result of several 
such simulations with varying initial L± are shown in Fig. 1111 
Each simulation used — 7r^/8 ~ 1.23 (so that Ux — U) 
as the interfacial free energy parameter, though system size 
varies as a function of L± . The parameters for this simulation 
were Ix 312, ly 104 (L 5.2687), and x/ 260. The 
final front speed was U = 0.00113 giving a scaling constant 
C = 0.17746. 
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FIG. 11. (color online) Measurement of the scaling constant 
C for the finger morphology shown in Fig. 1101 A proto-stripe 
morphology of wavelength L± is stabilized in a LBM simula- 
tion by adjusting the front speed U until a stationary profile 
is achieved. The fitted line gives an estimate of the minimum 
orthogonal stripe wavelength L™^^ for a given front speed U. 



C. Region of Stability 

We now have a theoretical prediction for the minimum 
and maximum orthogonal stripe wavelength as a function 
of front speed. These boundaries describe region of sta- 
bility for orthogonal stripes. Fronts moving faster than 
the speed at which the minimum and maximum stripes 
intersect JJ^^^ = 0.52 will not form orthogonal stripes. 
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(a)Initial configuration (b) Stable stripe 




(c)Merged via coarsening (d)Split via nucleation 



FIG. 12. Examples of simulations performed to determine 
the stability of orthogonal stripes of a given size formed by a 
front moving at a prescribed speed in order to map the region 
of stable stripe formation in Fig. 1131 The initial configura- 
tion (a) is similar to the one used to predict the lower bound 
of the stability region (see Fig. [TQ|) . and those observed in 
Figs. [1] m and [5] The simulation is run until the number of 
^-to-H-type interfaces intersecting the front and the outflow 
boundary changes, then the simulation is halted and classi- 
fied. If the number of interfaces at the front decreases, the 
simulation is classified as having merged (c) and the point in 
Fig. [13] gets a A symbol. If this number increases the sim- 
ulation is classified as split (d) and gets a V symbol. If the 
number of interfaces at the outflow boundary becomes 4 the 
simulation is classified stable (b) and the symbol is a circle 
with the radius proportional to the extremum value of the 
chemical potential at the front. 
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FIG. 13. Verification of minimum and maximum orthogo- 
nal stripe widths as a function of front speed for (j)in = 0. 
Filled circles are at points where simulations demonstrated 
stable stripe formation. Triangles are at points where stripe 
formation was unstable in the simulation. See Fig. [12] for a 
description of how the simulation results were obtained. The 
predicted stable region is above the minimum (dashed) and 
below the maximum (solid) lines, and corresponds well with 
the field of filled circles. See text for further discussion. 
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although orthogonal stripe may be formed by other in- 
fluences after the front passes: for example the presence 
of external walls, etc. as in [9]. 

We observe that the orthogonal stripes formed by a 
front in the simulation results shown so far fall in the 
region of stability, however we can test the bounds of the 
region more systematically. We do this by testing points 
in the vs. U parameter space for the formation of 
a stable stripe from the initial proto-stripe configuration 
shown in Fig. |12(a) If the point is inside the region of 
stability, the proto-stripe grows as shown in the time- 
lapse overlay image in Fig. 12(b) If the point is below or 



above the region, the proto-stripe will respectively merge 
or split; examples of which are shown in (c) and (d) of 
Fig. [12] If the point is to the right of the region, splitting 
and merging is followed by addition morphology changes 
until the stripes can reorient to become parallel to the 
front, however the simulation is classified and halted at 
the first morphology change. The final results of this se- 
ries of simulations is shown in Fig. [13] These simulations 
confirm the predicted region of stability for orthogonal 
stripes for symmetrically mixed (^in = 0) initial con- 
centration. Since our analytical predictions did not take 
into account the finite interface width we expect there 
will be deviation at small stripe widths and fast fronts. 
We see that the prediction for U^^^ = 0.52 is close to 
the maximum front speed capable of forming orthogonal 
stripes observed as U^^ ^ 0.3 in Fig. [l3]or U^"" ^ 0.2 
in Fig. [6] The predicted allowed stripes widths for slow 
moving fronts agrees very well with simulation results. 



However, a careful observer will notice that there is a 
slight discrepancy and our proto-stripe seem to allow for 
slightly smaller front speeds. We attribute this detail to 
the shape of the persistent stripe, which is different for 
both kinds of simulations (see Figs. [IQ] fc [T2]). The cur- 
vature in Fig. [To] will allow for a slightly faster coarsening 
of the stripe tip. 

We have now determined the characteristics of the 
ordered two-dimensional morphologies observed to be 
formed by phase separation fronts moving into mixtures 
of critical concentration: the orthogonal stripe morphol- 
ogy just presented, and the parallel stripe morphology-an 
essentially one-dimensional structure which we described 
and analyzed in previous work[14, 16]. The other ordered 
morphology we observed was an ordered droplet struc- 
ture which was not observed to occur for fronts moving 
into critical mixtures. 



VII. OUTLOOK 

In this paper we presented a survey of the morpholo- 
gies formed in the wake of a sharp phase-separation front. 
The resulting morphologies could be characterized as 
lamella in a parallel orientation with respect to the front, 
lamella in an orthogonal orientation and droplet arrays. 
We found that the selected morphology depended on the 
front speed, the volume fraction of the overtaken mate- 
rial, but also on the history of the system. If the front 
emerges from a homogeneous quench a depletion layer 
is typically formed and this will lead to the preferred 
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formation of lamella oriented parallel to the front. We 
saw, however, that sometimes defects will form and two 
systems under the same conditions-where the only dif- 
ference lies in the random initial conditions for the ho- 
mogeneous phase-separation-can instead form orthogo- 
nal lamella. By providing an unbiased initial condition 
we were able to determine the "preferred" , i.e. most sta- 
ble morphology. Using these results we were then able to 
present a state diagram as a function of the front speed 
and the volume fraction. We then examined in detail the 
formation of the orthogonal lamella for fronts moving into 
critical mixtures. We determined the range of allowed 
lamella sizes for a given front speed by analytically pre- 
dicting the minimum and maximum lamella which can 
be stably formed. This gave a prediction for the transi- 
tion between orthogonal and parallel lamella for critical 
mixtures which was within a factor of 3 of the observed 
transition point. 

The next step in this analysis will be to determine 
the allowed stripe widths for fronts moving into mix- 
tures with a minority and a majority phase in an ef- 
fort to predict the boundary of the "orthogonal stripes" 
region in Fig. [7]-now marked with an observed dashed 
line. This will require additional research, such as: De- 
termining how dynamical coarsening of stripe morpholo- 
gies is changed by having off-critical mixtures; this sub- 
ject is not well studied as dynamical scaling is typically 
studied in the context of homogeneous quench, which for 
off-critical mixtures results in droplet morphologies that 
coarsen due to Ostwald ripening [23]. Determining what 
effect off-critical mixtures ahead of the front will have on 
the nucleation of new stripe domains. Determining how 
the stripe morphology itself is altered by having an off- 
critical mixture. This last point may seem trivial at first 
glance, as conservation of the order parameter requires 
for stable stripes to form, the final width of the minority 
and majority stripes are a simple function of the mixed 
material concentration. We observe this to be true in 
our simulations for some distance away from the front, 
but this is not the case directly at the front where the 
dynamics of morphology selection occur. As shown most 
readily in Fig. SI but also elsewhere in this paper, there 
is a pinching off of the minority phase due to the pres- 
ence of a preferred contact angle induced by the control 



parameter front. These considerations will be published 
in future work. 

This paper raises a number of interesting issues. 
Firstly, can one predict a priori^ which morphology is 
most stable, and thereby provide analytical predictions 
for the state diagram? Secondly, what are the limits of 
metastability, i.e. what is the fastest orthogonal lamella 
structure that can be formed and what is the slowest par- 
allel lamella morphology that can be formed and how far 
can either encroach on the droplet states and vice versa? 
Thirdly, what is the region of stability for orthogonal 
lamella formed in off-critical mixtures, and by determin- 
ing this region can we predict the boundary between or- 
thogonal lemella droplet morphologies? Additionally, it 
would be very useful to more accurately determine where 
in the parameter space of this, and perhaps other simi- 
lar models, is the transition from parallel to orthogonal 
stripe morphologies. This transition has been noted sev- 
eral times, but a systematic study has not been done. 

So far we have been able to successfully analytically 
predict the size of parallel lamella structures. In future 
work we will refine our prediction for the extent of the 
orthogonal stripe region in the morphology diagram by 
considering the effect of the non-neutral wetting condi- 
tion at the front for off-critical volume fractions. An- 
other natural extension is to consider the system in three 
dimensions. In a future paper we will present an anal- 
ogous study which shows a slightly richer state diagram 
which includes cylinder arrays as well as three dimen- 
sional droplet lattices. Lastly this paper only considered 
diffusive dynamics. For many practical application it is 
important to include hydrodynamics effects, which can 
alter domain formation considerably. 
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